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Abstract 

We present the results of our investigation on the use of the two- 
body integrals to compute preliminary orbits by linking too short arcs 
of observations of celestial bodies. This work introduces a significant 
improvement with respect to the previous papers on the same subject 
13], 0. Here we find a univariate polynomial equation of degree 9 
in the radial distance p of the orbit at the mean epoch of one of 
the two arcs. This is obtained by a combination of the algebraic 
integrals of the two-body problem. Moreover, the elimination step, 
which in 13], 0 was done by resultant theory coupled with the discrete 
Fourier transform, is here obtained by elementary calculations. We 
also show some numerical tests to illustrate the performance of the 
new algorithm. 


1 Introduction 

This paper is related to the research started in i. a. where some hrst 
integrals of the two-body problem were used to write polynomial equations for 
the linkage problem with optical observations of asteroids and space debris. 
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In [3] the authors consider polynomial equations with total degree 48, that 
are consequences of the conservation of angular momentum and energy, and 
they propose two methods to search for all the solutions using algebraic 
elimination theory. The same equations were hrst introduced in |7] , but their 
algebraic character was not fully exploited at that time. In [3] the authors 
introduce for the same purpose new polynomial equations, with total degree 
20, using the angular momentum integral and a projection of Laplace-Lenz 
vector along a suitable direction. Both in [3] and in |1] the number of the 
considered equations is equal to the number of unknowns. 

Here we improve signihcantly the previous results by writing an overdeter¬ 
mined polynomial system (more equations than unknowns) which is proved 
to be generically consistent, i.e. the set of solutions in the complex held is 
not empty. For generic values of the data, by variable elimination, we obtain 
a system of two univariate polynomials of degree 10 with a greatest common 
divisor of degree 9. 

We also discuss the case where the oblateness of the Earth is relevant, so that 
we add the perturbation of the J 2 term to the Keplerian potential. In |2] 
and |3] this problem was faced by an iterative scheme, writing at each step 
polynomial equations with the same algebraic structure as in the unperturbed 
case (without J 2 effect). In this paper, the overdetermined system that can 
be written following the steps of the unperturbed case has the same algebraic 
structure but is generically inconsistent. However, we can use the iterative 
scheme mentioned above by neglecting one polynomial of the system. 

This paper is organized as follows. After recalling the linkage problem 
and some preliminaries on the two-body integrals (Sections [^[^, in Section]^ 
we discuss the polynomial equations that can be written for this problem, 
including the overdetermined system ([^ which is the object of this work. In 
Sectionwe show how the elimination steps can be carried out, and we prove 
the consistency of system (|^. In Section]^ we discuss the spurious solutions 
of and illustrate two methods to discard them. Section is devoted to 
the linkage problem with the J 2 effect. Some numerical tests are presented 
in Section In the Appendix we discuss a simple way to hlter the pairs of 
attributables to be linked with this algorithm. 
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2 Linkage of too short arcs 


We consider objects moving in a central force field. Let ns fix an inertial 
reference frame, with the origin at the center of attraction O, which is the 
center of the Snn (Earth) in the asteroid (space debris) case. Assnme the 
position q and velocity q of the observer are known fnnctions of time. We 
describe the position of the observed body as the snm 

r = q + pe^, 


with p the topocentric distance and the line of sight nnit vector. We 
choose spherical coordinates {a,S,p) G [—vr, tt) x (—7r/2,7r/2) x M+, so that 

= (cos 5 cos a, cos 5 sin a, sin 5). 

A typical choice for a, S is right ascension and declination. The velocity 
vector is 


f = q + pe^ + p{a cos5e" + ^e"^), p, d, 5 G M, p G 


where p, pdcos5, p6 are the components of the velocity, relative to the 
observer, in the (positively oriented) orthonormal basis {e^, e", e^}, with 


e 


a 


(cos 5) ^ 


deP 
da ’ 


e^ = 


deP 


Let {ti,ai,6i) with z = 1... m, m > 2, be a short arc of optical obser¬ 
vations of a moving body, made from the same station. If m > 3, we can 
compnte a, 6, d, 6, d, 6 at the mean time t = A U by a qnadratic 
£t. From these qnantities we can try to compnte a preliminary orbit. When 
the second derivatives are not reliable dne to errors in the observations (or 
not available, if m = 2) we speak of a too short arc (TSA) and, to compnte 
a preliminary orbit, we have to add information coming from other arcs of 
observations. This is a typical identihcation problem, see [6]. 

In any case, it is possible to compnte an attributable 


A = {a, 6, d, (5) G [— TT, tt) x (—7r/2,7r/2) x 


representing the angnlar position and velocity of the body at epoch t (see 
0.0). The radial distance and velocity p, p are completely nndetermined 
and are the missing qnantities to define an orbit for the body. 

In this paper we deal with the linkage problem^ that is to join together two 
TSAs of observations to form an orbit htting all the data. 
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3 First integrals of Kepler’s motion 


We consider the first integrals of the equation of Kepler’s problem 


r = 


as functions of the unknowns p, p. The angular momentum is the polynomial 
vector 

c(p, p) = r X f = Dp + Ep^ + Fp + G, 

with 


D = q X e^, E = e'’ X e-*-, F = q x e-*- + x q, G = q x q, 


where we have set 

e"*" = d cos 5e“ + 

The expression of the energy is 

- A, (1) 

2 |r| 

where 

|r| = (p2 + |q^ + 2pq■e0'/^ ( 2 ) 

|f|^ = p^ + |e-^|^p^ + 2q ■ e^p + 2q ■ e-^p+|qp. (3) 

The Laplace-Lenz vector L is given by 

pL(p,p) = f X c — p— = ^|f|^ — (4) 

with |r|, |i-|2 as in g, and 

f-r = pp + q ■ e^p + (q ■ + q ■ e-‘“)p + q ■ q. 

Moreover, the following relations hold for all p, p: 

c ■ L = 0, p^|L|^ = + 2T|c|^. (5) 


Expressions 0.0 are algebraic, but not polynomial, in p, p. However, we 
can introduce a new variable u G M, together with the relation \y\^v? = p^ 
and we obtain 

8 = ~ ~ ( 1 ^ 1 ^ ~ ~ 
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that are polynomials in p, p, u. 

For later reference we also introduce the quantity 

K = ~ (®) 

4 Polynomial equations for the linkage 

We use the notation above, with index 1 or 2 referring to the epoch. Let 

jXj (oj, 6j , Q.j 5 ), j 1)2 

be two attributables at epochs ij. We consider the polynomial system 

Ci=C 2 , Li = L2, 81=82, Mi|rip = p^ M^|r 2 p = p^ ( 7 ) 

in the 6 unknowns 

{Pi, P2, Pi, P2 ,Ui,U2)- 

System Q is dehned by the vector of parameters 

(A,^2,qi,q2,qi,q2), 

and we shall discuss properties which hold for generic values of them. More¬ 
over, Q is composed by 9 equations with 6 unknowns. However, due to 
relations (|^, 2 equations can be considered as consequences of the others, so 
that we are left with a system of 7 equations with 6 unknowns. 

Assume the attributables Mi,M 2 refer to the same observed body, the 
two-body dynamics is perfectly respected, and there are no observing er¬ 
rors. Then the set of solutions of Q is not empty. Taking into account 
the observational errors, and the fact that the two-body motion is only an 
approximation, system Q turns out to be generically inconsistent, i.e. it has 
no solution in C. 

We search for a polynomial system, consequence of Q, which is generi¬ 
cally consistent, with a hnite number of solutions in C, and which leads by 
elimination to a univariate polynomial equation of the lowest degree possible. 

Introducing relations u‘j\rj\‘^ = (j = 1,2) for the auxiliary variables ui, 
U 2 corresponds to the squaring operations, used in B. 01 to bring the selected 
algebraic system in the variables (pi, p2, Pi, P 2 ) into a polynomial form. Since 
these operations are responsible of the high total degree of the resulting 
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polynomial systems (48 and 20 respectively), in writing the new eqnations 
we try to cancel the dependence on both mi, U 2 by algebraic manipnlations 
of the conservation laws. First we shall consider the intermediate system 

Cl = C2, p(Li - La) = (£^i - £^2)r2, M?|rip = (8) 

where U 2 does not appear, which is still inconsistent; then we shall take into 
acconnt the system 

Cl = C2, (Ki - Ka) X (ri - ra) = 0, (9) 

where also Ui does not appear, whose consistency is proven in the next sec¬ 
tion. 


5 Elimination of variables 

In this section we show that generically system (|^ is consistent. In particnlar, 
by elimination of variables, we shall end np with two nnivariate polynomials 
of degree 10 in the range pa, whose greatest common divisor generically has 
degree 9. A similar procednre can be carried ont by eliminating all the 
variables bnt pi. 


5.1 Angular momentum equations 

The conservation of angnlar momentnm gives ns 3 polynomial eqnations that 
are linear in pi,p 2 , and qnadratic in pi, p 2 - Therefore, it is natnral to nse 
these eqnations to eliminate the radial velocities, as done in [3], [Ij. These 
eqnations can be written as 


Dipi — D2P2 — J (pi, P2), 


( 10 ) 


with J a vector whose components are qnadratic polynomials in pi, pa. Fol¬ 


lowing [3] we project (10) onto the vectors Da x (Di x Da) and Di x (Di x Da) 
and obtain pi, pa as qnadratic polynomials in pi, pa. With these expressions 
of pi. Pa we have 

Ac X (Di X Da) = 0, (11) 


with Ac = Cl — Ca, whatever the valnes of pi, pa. The projection of (10) onto 
Di X Da allows ns to eliminate the variables pi,p 2 and yields 


q{Pi, P 2 ) — Q 2 fiPi + qiflPi + qo, 2 P 2 + Qo,iP 2 + <?o,o, 


( 12 ) 
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where the coefficients qij depend only on the attributables and on the posi¬ 
tion and velocity of the observer at epochs ti, t 2 . 

In the following we shall consider the quantities introduced in Section as 
function of pi, p 2 only, by the elimination of pi, p 2 just recalled. 

5.2 Bivariate equations for the linkage 

By subtracting £ 2'^2 to both members of the Laplace-Lenz equation, and 
using the conservation of energy, we obtain 


P-Li — S1V2 — p^2 ~ £-2^2- 


Equation (13) can be written 


= 0 , 


where we have set 


(13) 


(14) 


Ak = Ki - K 


2, 


= ri - r2. 


with K as in (§. Note that the variable U 2 does not appear in did] ). 
We can also eliminate Ui by cross product with A^: 


AkX Ar = 0 . (15) 

For brevity we set ^ = A^ x A,., and we note that 

^ - IriHri X r 2 - (fi ■ ri)fi x A^ ( 1-2 ■ r 2 )f 2 x A^. (16) 

Remark 1 . By developing the expressions of ri x r2, fi x A^, f2 x A^ as 
polynomials in pi,p 2 we obtain that the monomials in ^ with the highest total 
degree, which is 6, are all multiplied by ej* x e^. 

In the following section we shall prove that for generic values of the data the 
system 

q = 0, ^ = 0 (17) 

is consistent, that is the set of its roots in C is not empty. 

Note that, if g = 0, the vector ^ is parallel to the common value Ci = C 2 of 
the angular momentum. 
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5.3 Consistency of the equations 

The proof of the consistency relies on some geometrical considerations. In 
particular it is relevant to check whether the angular momentum vector is 
orthogonal to the line of sight. We introduce the quantities 

Cij=Ci-eP, z,j = l,2. 


More explicitly we have 

Cii = qi X + qi x ■ qi, 

ci2 = X + qi x ■ e^pi(pi, ps) + 

+ (e^ X qi + qi X e^) ■ e^pi + qi x qi ■ e^, 

and similar expressions for C 22 , C 21 . In particular, equations cn = 0 and 
C 22 = 0 represent straight lines in the plane pip 2 , while C 12 = 0 and C 21 = 0 
give conic sections, see Figure [T] 



Figure 1: In the plane pip 2 , for a test case, we draw the curves g = 0 (black), 
C 12 = 0 and C 21 = 0 (light gray), and the straight lines Cn = 0, C 22 = 0 
(dashed), and pi = Pi,p 2 = p '2 (dotted). 









Consider the point C = {p'l, p'i) defined by Cn = C 22 = 0, so that 


Pi 


qi X qi ■ 

X . qi ’ 


// 


P2 


q2 X q2 • 

62 X ■ q2 ■ 


(18) 


In Lemma we shall prove that C lies on the conic g = 0 and is the only 
point where both angular momenta Ci, C 2 vanish. The straight line p 2 = P 2 
generically meets g = 0 in another point P 2 = {PiiP’i)^ where the angular 
momenta do not vanish. Similarly, the straight line pi = p'l generically meets 
g = 0 in another point Pi = (^ 1 ,^ 2 )! where the angular momenta are not 
zero, see Figure[T] For g = 0, the vector ri x r 2 gives the direction of Ci = C 2 , 
therefore from the equations ri x r 2 ■ = 0, j = 1, 2 we obtain 


p'l 


qi X q 2 ■ 

e? X ■ q 2 ’ 


P2 


qi X q 2 ■ 

e? X ■ qi ■ 


(19) 


Note that 


g 2 ,o = -El • Di X D 2 = -(e^* X ■ qi)(e(’ X ■ qa). 


so that g 2 ,o 7 ^ 0 implies that both and p'{ are well defined. 

In a similar way we obtain 

go, 2 = E2 ■ Di X Da = (e^ X ■ q2)(e(’ X ■ qi), 

so that go ,2 7 ^ 0 implies that both p '2 and P 2 are well defined. 
Generically we have 

^2,0, go,2 7^ 0. 


( 20 ) 


Lemma 1. If (20) holds, then the point C = {p'l, P 2 ) given by Cn = C22 = 0 
satisfies q{Pi,P2) = 0 . Moreover, in C we have Ci = Ca = 0 and C is the 
unique point in the plane pip 2 where both angular momenta vanish. 

Proof. Using relation r, ■ D,- = 0 we obtain 


Cj X D,- = -(fj • D,)r,-, 


J = l,2. 


Moreover, condition (20) yields e? X q, ^ 0, so that Yj 7 ^ 0. 


From relations 


G ■ Dj = 


-c ■ • e^ 


9 








we have 


( 21 ) 


Cj X Dj = 0 if and only if Cjj = 0 
for j = 1, 2. Therefore cn = C 22 = 0 implies 


Ac ■ Di X D 2 = 0, 


that together with ([II| gives 


Cl = C 2 . 


Finally, relations (21), (22) imply Cj = 0, j 
ately follows from the definition of C. 


( 22 ) 

1,2. The nniqneness immedi- 

□ 

□ 


Lemma 2. In the point C = {p'i, p'i) generically we have ^ 7 ^ 0. 
Proof. By Lemma if (20) holds, we have Ci = C 2 = 0 in C, so that 

p(Li — L2) — (Ti — £2)^2 = p(-, —r — -j—r) — (£1 — £2)^2 


P2 




^ , f\- \2 I• |2\ 

= -p -—r + -(|r2| - |ri| )r2. 


Til 


Therefore we have 


kx X Ac = -(|fi| 


|r2| jri X r2. 


(23) 


We show that the right-hand side of (23) does not vanish in C. In fact, by 
projecting ri x r 2 onto qi, q 2 we obtain 


qi ■ ri X ra = p" [p'iiei X ■ qi) - qi x qa ■ e{] , 
qa ■ ri X ra = Pa [p'li^i x ■ qa) - qi x qa ■ e^] , 


and the expressions in the brackets vanish only if 


Pi ~ Pi 5 P 2 ~ £ 2 - 

Moreover, p" = p 2 = 0 occurs only if q^ x q^ ■ = 0, j = 1,2. Thus 

ri X ra generically does not vanish. Using Ci = 0 projected e.g. onto e;;*-, we 
can prove that in C the quantity |fip depends only on the data Tli,qi,qi 
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at epoch ii. In fact through this relation we can hnd an expression for pi 
depending only on ^i,qi,qi. A similar result holds for |f 2 p. We observe 
that the |fjp do not vanish individually. Indeed, if relations < X ■ 7^ 0 ; 

j = 1,2 hold, then fj 

P does not vanish. We conclude that the 


— A nP 


Lj = pjGj + PjBj + q^ 7 ^ 0 for any choice of Pj,pj- 
Therefore generically also |fip — |f 2 


point C is not a solution of (17). 


□ 


□ 


Lemma 3. Assume q = 0. Then ^ = 0 is generically equivalent to 


^e? = 0 
c ■ e? 7 ^ 0 


or 


^ = 0 

c ■ 7 ^ 0 


(24) 


Proof. Assume 


does not hold. Clearly relations ^ 


= ^ 


= 0 


are necessary to have ^ = 0. Then we have c ■ = c ■ = 0. If (20) 


holds. Lemma [T] implies Ci = C 2 = 0 and, by Lemma generically we have 
^ ^ 0 . Viceversa, assuming ^ 7 ^ 0 we obtain that each system in ( [24) ) is 
incompatible, because for g = 0 we have ^ x c = 0. 


□ 


□ 


Lemma 4. Assume relation 


A, ■ e? X 


7^0 


(25) 


holds. Then ^ = 0 is equivalent to 

^ = 0 and ^-62 = 0 . 


Proof. Relation ^ = 0 holds trivially. Moreover, since (25) is satisfied, 

the vectors e^’, e^, A^ are linearly independent. 


□ 


□ 


We introduce the polynomials 

Pi = ^ ■ e?, P2 = ^ ■ e^. 
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(26) 


By Remark both pi and p 2 have total degree 5 in the variables pi, p 2 - 
We consider the system 

g = Pi = P2 = 0. 

We are now ready to state the main resnlt. 


Theorem 1. Generically, system (11) is consistent and can be reduced to a 
system of two univariate polynomials Ui,U 2 whose greatest common divisor 
has degree 9. 


Proof. By Lemma we only need to prove consistency of system (26). First 
we show that system ( |I^ has at least 9 solntions. In fact Lemma [^implies 
that both systems q = pi = 0 and q = p 2 = 0 dehne generically 10 points. 
Moreover, for g = 0 relation Cn ^ 0 discards the points Pi, C, while relation 
C 22 7 ^ 0 discards the points P 2 , C. By Lemma we know that generically C 
is not a solntion of ^ = 0 , hence it does not solve either pi = 0 or p 2 = 0. 
On the other hand, we have pi(Pi) = P 2 (P 2 ) = 0. We show that pi(Pi) = 0. 
If ^(Pi) = 0 the resnlts trivially holds. If ^(Pi) 7 ^ 0, by Lemma we have 
c(Pi) 7 ^ 0, so that c(Pi) II ^(Pi), and from cii(Pi) = 0 we obtain pi(Pi) = 0. 
In a similar way we can prove that P 2 (P 2 ) = 0. Therefore we are left with 9 
solntions for each system. We show that they are the same ones. In fact, by 
Lemma 1 ^ solntions of dl?) ) mnst satisfy pi = P 2 = 0. Moreover, generically 
we have Pi(P 2 ) 7 ^ 0 and P 2 {Pi) 7 ^ 0, so that both Pi and P 2 are not solntions. 
Therefore we have exactly 9 solntions. 

Let us consider the univariate polynomials 


u, = res(p,,(/,p,), 1 = 1,2, 


(27) 


that are the resultant of the pairs pj, q with respect to pi (see pQ). 

The root p 2 = p '2 of Ui and the root p 2 = P 2 of U 2 must be discarded because 
they correspond to the points Pi, P 2 for the polynomials pi,p 2 respectively. 
We consider 

Ui _ U2 


Ui = 


P2 - P2 ’ 


U2 = 


P2- P2 


By the discussion above, these polynomials have degree 9 and must have the 
same roots: in particular, up to a constant factor, they both correspond to 
the greatest common divisor of Ui and U 2 . 


□ 


□ 
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5.4 The univariate polynomials Ui, U 2 


We explicitly perform the elimination step to pass from system q = pi = 
P 2 = 0 to 

Ui = U 2 = 0, 


where Ui,U 2 are the two univariate polynomials introduced in (27). For this 
purpose we produce an equivalent system q = pi = p 2 = 0 where the pj are 
linear in one variable, say pi. 

Assume q is not degenerate, i.e. ^ 2 , 0 , qop 7 ^ 0. Indeed, here we use ^ 2,0 7 ^ 0 
only, while go ,2 7 ^ 0 is necessary for the similar construction relative to p 2 - 
We write 


q{pi,P2) = ^iqp\p2 

i,j=0 


2 


= > Mp2)Pi, 


h=0 


where 


^ 0 (^ 2 ) — qoppl + qo,iP2 + qo,o, bi — qifl, 62 — q2,o- 

Moreover we write 

5 4 

Pi{pi,p2) = (28) 

i,j=0 h=0 

5 5 

P2{P1,P2) = ^Pf]p\p’2 = ^(^2^{p2)p'l. (29) 

i,j=0 h=0 


for some polynomials au^h whose degrees are described by the small circles 
used to construct Newton’s polygons of pi,P 2 in Figure From g = 0 we 
obtain 


pJ = ft,Pi + 7ft, ft = 2,3,4,5 


(30) 


where 


and 


ft = 4 . 

02 


72 = 


4+1 — 44 + 7/u 7/ft+i — 472, h — 2,3, 4. 


Inserting (30) into (^28l, (29) we obtain 


Pj{pi, P 2 ) = ai,i(p 2 )pi + ai,o(p 2 ), j = 1, 2 


(31) 
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where 


4 

hi,! = ai,l + ttl,hPh, 

h=2 

5 

02,1 = 02,1 + 0,2,h(3h-, 

h=2 



0 2 4 6 

Pl 



0 2 4 6 

Pl 


4 

Ol,0 = Oi,o + 0‘l,hlh-, 

h=2 

5 

02,0 = 02,0 + 0.2,hlh- 

h=2 



Pi 



0 2 4 6 

Pi 


Figure 2: Newton’s polygons Q, Pj, Pj for the polynomials g, Pj,pj, j = 
1,2. The nodes with circles correspond to (multi-index) exponents of the 
monomials in pj, pj] the nodes with asterisks correspond to exponents of the 
monomials in q. 


In Figure bottom left and right, we draw Newton’s polygons of pi, p 2 , 
which also describe the degrees of the polynomials dk,h- 


From (31) we get two expressions for pi. 


Pl = 


Ol,0 

di,i 


Pl = 


02,0 

02,1 
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By substituting these expressions into q we obtain two univariate polynomials 
of degree 10 in the variable p2- 


— q2flOb\fl — 

^2 = '?2,0^2,0 ~ <?1,002,002,1 + ^oO^y- 

Using the properties of resultants we hnd that 

^1 = ^2 = q %^ 2 - 


5.5 Non-degeneracy conditions 

We list below the conditions on the data .Ai, .4,2, qi, q2, qi, q2 that we used 
in the previous sections. These conditions generically hold. 


1 . El ■ Di X D2,E2 ■ Di X D2 7^ 0 , so that g is a quadratic polynomial 
both in Pi and p2- An interpretation of these relations is given in [ 3 ]. 
Moreover, these conditions imply: 


i) Di X D2 7^ 0 , so that we can compute pj = Pj{pi, P2) from system 
( 10 ). This condition also implies Di,D2 7^ 0 , which in turn yield 
ri,r2 7^ 0 for all pi, P2] 


ii) ei X e;j*“ ■ qi, 62 X 62 ■ q2 7^ 0 , which are used to dehne p'i, P2 
respectively; 

hi) 61 X 63 ■ qi, 61 X e2 ■ q 2 7^ 0 , which are used to dehne P2, p) respec¬ 
tively. 


2 . qi X q2 ^ 0 , P'i ^ Pp P'l ^ ~ These conditions imply 

ri X r2 7^ 0 , and also ri, r2, 7^ 0 , for all pi, P2- 


3 . 

4 . 


X 


■ qj 7^ 0, so that rj 7^ 0 for all pj, pj, with j = 1, 2. 


Aq-ei X 63 7^ 0 , so that {ei, e^, A^} forms a basis of M^. This condition 
implies x 7^ 0, so that the maximal total degree for the components 
of equation ( 15 ) is 6. 


5 . Pi(P2),P2(-Pi) 7^ 0 , so that Pi, P2 are not solutions of ([TT]). 
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6 Selecting the solutions 


After computing all the solutions of © we can select the ones with both 
entries pi,p 2 real and positive and compute the corresponding values of 
Pi, P 2 - However, since equations ([^ impose only some of the laws of the 
two-body dynamics, we expect that some of the remaining solutions yield 
vectors (pi,P 2 ,Pi,P 2 ,^ 2 ) which do not solve ([^. 

In this section we characterize the spurious solutions, and propose some 
algorithms to select the good ones. To decide which solutions can be accepted 
we suggest to use one of the two methods introduced in 0 . HI- They take 
into account the errors in the observations, which can be represented by 
4x4 covariance matrices Ti, r 2 of the attributables. We recall that the hrst 
method relies on the computation of a norm referring to some compatibility 
conditions, see |3], Section 5. The second method requires to compute a 
covariance matrix for the candidate preliminary orbits, which is used for 
the attribution algorithm, see [1], Sections 7, 8. The new formulas for the 


covariance matrix are provided in Section 6.2 


6.1 Spurious solutions 

We consider real solutions of ([^ that do not solve the intermediate system 
(|^, and solutions of Q that do not solve Q. In the hrst case the spurious 
solutions do not satisfy 


A.K ■ Ar -I- 



= 0 . 


Concerning the second case we prove the following: 

Proposition 1. Each real solution of ^ which does not solve 0 fulfills 


Li — L 2 I — 2. 


Proof. Using Ci = C 2 , the second equation in can be written as 


(Li — L 2 ) — p 


(|Lil 

2 

— 

IL 2 I 

P) 

2|c| 

2 


r2, 


(32) 
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r2 T2 




Figure 3: Left: case |Li| > IL 2 I. Right: case |Li| < IL 2 


where c is the common value of the angular momentum. The conservation 
of Laplace-Lenz vector and energy is in general not guaranteed. However, 


we note that if Ci = C 2 and |Li| = IL 2 I, then ([^, (32) imply Li = L 2 and 
= ^ 2 . 

If |Li| > IL 2 I, the vectors Li — L 2 , r 2 have the same orientation. Passing to 
the norms in (32) and substituting 



c 

2 

/i(l + 

L2 

COS 6*2) ’ 


T2I = 

where 62 is the angle between L 2 and r 2 , we obtain 

Li — L2 = 



|Li| 

2 _ 

IL 2 I 

2 

2(1 + 

IL 2 I 

cos 

02 ) 


Using relation 


|Li|2-|L2p = IL 1 -L 2 


2 |Li — L 2 IIL 2 I cos 62 


and rearranging the terms we get |Li — L 2 I = 2. 

If |Li| < IL 2 I the vectors Li — L 2 , r 2 have opposite orientation, and we 
obtain |Li — L 2 I = —2, which is impossible. In fact in this case we have (see 
Figure]^ 


— IL 2 P — |Li — L 2 P T 2|L]^ — L 2 IIL 2 I cos(7r — O 2 


□ 

□ 

Since system ([^ is generically inconsistent, the spurious solutions discussed 
in Proposition^ usually do not occur. 
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6.2 Covariance of the solutions 

We introduce the vectors 


A — (Ai, A. 2 ), 


It — (Pn Pi, P2, P2) 


and the map 


$(R,A) 


Cl - C2 


Assume we have a pair of attributables A = (Ai, A. 2 ), with covariance ma¬ 
trices Fi, r 2 , at epochs ti, ^ 2 - For each vector R such that 


$(r,A) = o, 


det 




(R,A)^0 


there exists a map A 1 —)■ R(A) with $(R(A), A) = 0, dehned in a neighbor¬ 
hood of A, with R(A) = R. Following |1], we consider the map defined 
by $ = o Tau^, where T^tT is the transformation from attributable coor¬ 
dinates Eatt to Cartesian coordinates Ecar at the two epochs. Then we can 
use the same scheme as in |1] (Section 7) to compute the covariance matrix 
of the Cartesian coordinates at epoch ti through the formula 


f(i) 

car 


dA 


r, 



T 

r A — 

r Fi 

0 

y dA \ 

•) 

A 

0 

F2 J 


with the derivatives evaluated at A = A. 

The only differences with respect to IH are in the term which in this 

6^ Onjcar ' 

case is given by 


dErn, 


where we use the hat map 


9 (mi,M2,%) 


-ri 

ri 

r2 

-r2 

dpi 

dpi 

dpi 

dpi 

. dri 

dri 

dr2 

1 



0 

-Ms 

U2 

= U 1 —u = 

U 3 

0 

-Ml 


-U 2 

Ml 

0 
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To compute the derivatives of pi = x we take advantage of the 

following relation, which holds at the solutions of pi = 0: 

dpi 1 dpi 

^Ecar Pi dE^ar 

Here pi = A.k x A^ ■ vi, with vi = ri — qi, and 
dn* 1 

^ = -(|f2|^ - |fiHqi X ra - (fi X A^ ■ vi)fi 

OVI Z 

-(ri ■ ri)[-qi x fi + ra x fi] + (fa ■ ra)[-qi x fa + ra x fa], 

^ = -(ri X ra ■ vi)fi - (fi X Ar ■ vi)ri - (fi ■ ri)A^ X vi, 

OTi 

dn* 1 

-pp = -;r(|f2p - li-inqi X ri - (fi ■ ri)fi X Vi 

OV 2 z 

+ (fa X A^ ■ vi)fa + (fa ■ r 2 )fa x Vi, 


dpi 

dr2 


(ri X ra ■ Vi)fa + (fa X A^ ■ Vi)ra + (fa ■ ra)A^ x Vi. 


7 Including the J 2 effect 


If we want to compute preliminary orbits of space debris, we have to take 
into account the effect of the oblateness of the Earth (Ja effect), which gives 
rise to a perturbation of the two-body dynamics. Following [2] and [1] we 
can use an iterative scheme to include the Ja effect in the determination of 
preliminary orbits with the two-body integrals. We consider the equations 

= Ca, i^iLi = La, 81 = 82, u 1 \ri\'^ = p, ul\r2\‘^ = 

where the rotation matrices 

Rc = -Rao) Rl = RJi+AujRAnR'^uii (33) 

are dehned through the angles 

Aca = UJ2 — All = fla — fli. 


In (33) we use Rg to denote the rotation by the angle 9 around the axis 
dehned by the vector v. 
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Following the same steps of Section we consider the intermediate equation 

— 8x^2 = / iL2 — 82X2-, 


which can be written 

^A- + (^|rir-Wi)^r =0, 

where 

= RlKi - Ka, = RlVi - ra. 

Then we can eliminate the dependence on ui by vector product with A^. 
We end up with the system 


RfX^x Ca, 


A^ 


X A^ = 0. 


(34) 


Note that system (34) is not polynomial due to the presence of the rotation 
matrices Rc, Rl, that depend on the orbital elements. We can search for 
solutions of (34) by considering the solutions of ([^, i.e. assuming Rl = 
Rc = I, as hrst guesses of the iterative method. Inserting these solutions 


into Rl and Rc, system (34) becomes polynomial, and we can solve it like in 


the unperturbed case. We iterate the procedure and consider the solutions 
at convergence, if any. Some care must be taken in selecting the solutions at 
each iteration. 


Similarly to (16) we have 


Ar X A^ 




|fi|^)i?Lri X ra - (fi ■ Ti)RLri x A^ + (fa ■ ra)f2 x A 


By developing the expressions of Rl^i x ra, Rl^i x A^ 
nomials in pi,pa we hnd that the monomials in Ar x 


^, fa X A^ as poly- 
A^ with the highest 


degree (i.e. 6) are all multiplied by Rl^I x e^. Then we can project the 
second equation in ( [^ onto the vectors Rl^i, ea, to obtain two polynomial 
equations Pj{pi, P 2 ) = 0, j = 1,2. As before, pi and pa have degree 5. The 
range rates can be eliminated from equation RcC\ = Ca, to get a quadratic 


polynomial q analogous to (12) 


An important difference with respect to Section is that in this case the 
system pi = pa = ? = 0 is generically inconsistent, as can be checked by a 
numerical test. Nevertheless we can choose either pi = g = 0 or pa = g = 0 
as polynomial equations for the linkage. 
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Pi Pi 


Figure 4: Left: Intersection of curves defined hj q = pi = p 2 = 0 for (101955) 
Bennu. Right: the same for (100000) Astronautica. 

8 Numerical tests 

We describe the results of two numerical tests, one for a near-Earth asteroid, 
the other for a main belt. The hrst object is (101955) Bennu. We link 5 
observations made on September 11, 1999 together with 11 observations made 
on March 30, 2000. After discarding non-real and non-positive solutions, we 
are left with only one pair 

(pi,p2) = (0.04379,0.27132). 

The second object is the main belt asteroid (100000) Astronautica. We link 
10 observations made on October 17 and 19, 2003 together with 4 observa¬ 
tions made on May 9, 2005. In this case we are left with the pairs 

(pi,p2) = (0.05240,0.02179), (0.83897,1.08648), (1.09052,3.04874). 

The third pair is discarded because it yields an unbounded orbit at epoch 0- 
The values of the penalty (see 0) for the hrst and second solution, computed 
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asteroid 

MJD 

a 

e 

/ 

n 

u 

i 

(101955) 

51432.4 

51633.4 

1.1316 

1.1306 

0.2058 

0.2039 

5.9895 

5.9895 

2.1229 

2.1229 

65.2425 

66.4329 

303.8601 

107.6598 

( 100000)1 

52930.2 

53499.3 

1.0365 

1.0367 

0.0277 

0.0305 

1.0298 

1.0298 

202.3124 

202.3124 

275.0221 

308.8393 

268.4023 

73.7641 

( 100000)2 

52930.2 

53499.3 

1.8725 

1.8754 

0.0768 

0.0862 

20.8665 

20.8665 

186.7970 

186.7970 

197.8094 

197.7172 

344.8383 

198.9997 


Table 1: Keplerian elements found for (101955) Bennu, and for (100000) 
Astronautica. 


by attribution, are X 4 = 1224012.479 and X 4 = 1-497 respectively, therefore 
we select the second solution. 

In Figure]^ we draw the curves g = 0 (black), pi = 0 (light gray), P 2 = 0 (dark 
gray). The dashed straight lines correspond to cn = 0 (vertical), C 22 = 0 
(horizontal). The computed pairs (pi,p 2 ) are marked with asterisks. 

In Table [T] we show for both asteroids the Keplerian elements at the two 
mean epochs ti,t 2 that we find with this method. For asteroid (100000), we 
indicate with the labels 1, 2 the two bounded orbits that we can compute. 
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A Appendix 

We suggest a way to discard pairs of attributables which are not likely to 
fulhll the requirements to be linked. This filter is based on the angular 
momentum equation Ci = C 2 . 

A.l Filtering pairs of attributables 

The equation g = 0 represents a conic Q in the plane pip 2 - We can decide 
to accept only pairs (pi,p 2 ) inside a square TZ = [pmin, Pmax] x [pmin, Pmax], 
with 0 < pmin < Pmax- We clieck whether the conic Q does intersect 71. 
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First consider the case Q is unbounded (hyperbola or parabola). In this case 
it is sufficient to check whether the conic intersects the boundary dTZ. 

If Q is bounded (ellipse or circle) this check is not enough; it can happen 
that Q n 7^ 7 ^ 0, with Q lying totally inside TZ. 

Here we sketch an algorithm for this hlter. Consider the four straight lines 

C {{PIt P2) ■ Pj Pmiri\ 1 ^j +2 {{.Plj P2) ■ Pj Pmax\ 1 fo^ j 1)2. 

First compute the intersections of the conic with each line j = 1... 4, if 
any. If an intersection lies on the segment of some belonging to dTZ we 
accept the pair of attributables and continue the linkage procedure. If not, 
we check whether Q is bounded or unbounded. 

If Q is unbounded we discard the pairs of attributables. If Q is bounded, we 
can check whether Q is totally inside TZ by computing the coordinates of the 
center of the conic. If the center is inside TZ we accept the pair provided Q 
has no intersections with any r^, otherwise we reject it. 
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